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GLOSSARY 


A 

B 

ci 

c(z) 

D,  E.  F,  G,  B(K.z) 
f 

g(x,Zl£,b>) 

h(x) 

K 


inversion  amplitude  (10) 
inversion  amplitude  ratio  (28) 
velocity  below  reflector  (35) 
reference  velocity  (1) 
various  integrals!  see  Appendix  A 
frequency  (53) 

Green's  function  for  c(z)  medium  (2) 
cylindrical  surface  (29) 
ray  parameter  (3) 


k, (E,z) 

n(z)  =  c(0)/c(z), 

R 

Rj 

S 

US<5;to>) 

v(x.z) 

w 

x  *  (x.y) 
z 

u(x» z) 

0 

Y 

T» 


Jn*(z)-K*  (6) 

the  index  of  refraction  (7) 
reflection  coefficient  (34) 
reflection  coefficient  (26) 
abbreviation*  see  (33) 
observed  scattered  field  (2) 
observed  scattered  field  (55) 
velocity  (1) 

the  inversion  operator  (10) 
horizontal  cartesian  (1) 
vertical  cartesian  (1) 
unknown  perturbation  in  velocity  (1) 
reflectivity  function  (26) 
abbreviation*  see  (32) 
abbreviation*  see  (34) 
bandlimited  delta  function  (26) 


i 


jump  in  e  (35) 

cartesian*  for  observation  point  at  z  *  0  (2) 

offset  (4) 

offset  (11) 

travel  time  (8) 

phase  function  (14) 

circular  frequency  (2) 


ii 


I 


0  r  i 


i 


Carter  and  Frazer  [1984],  and  Bleisteln  and  Gray  [1984]  (henceforth  BG) , 

I 

present  inversion  algorithms  which  include  the  effect  of  a  stratified 
reference  velocity.  c(z).  Those  papers  did  not  address  the  question  of 
obtaining  accurate  values  of  the  reflection  coefficient!  this  is  the  issue 

I 

treated  here.  Thus,  in  the  language  of  Bleistein,  Cohen  and  Hagin  [1984], 
(henceforth  BCH) ,  the  earlier  algorithms  provided  structural  inversions,  the 
^  location  of  the  sub-surface  layers*  whereas  the  present  algorithm  also 

provides  an  accurate  estimate  of  the  reflectivity  function,  which  depicts 
the  reflectors  and  provides  an  estimate  of  the  reflection  strengths  across 
the  layers. 

Since  we  employ  a  perturbation  assumption  (the  "Born  Approximation"),  the 
constant  reference  speed  inversion  first  described  in  Bleistein  and  Cohen 
[1979a]  and  reviewed  in  BCH,  is  often  not  adequate  at  depth.  Although 
recursive  use  of  the  algorithm  is  possible  and  although  the  results  can  be 
significantly  enhanced  by  suitable  pre-  or  post-processing  (e.g.,  see  Hagin 

I 

and  Cohen  [1984]),  extension  of  the  perturbation  method  to  a  stratified 
reference  profile  is  highly  significant.  It  is  far  more  likely  that  the 

actual  velocity  function  can  be  well  approximated  by  a  stratified  reference 

I 

velocity  than  by  a  constant  one,  which  in  turn  enhances  the  validity  of  the 
perturbation  assumption  and  the  inversion  results.  See  BG  for  further 

discussion  of  this  point. 

I 


The  algorithm  presented  here  has  the  same  structure  as  the  BG  algorithm 
and  hence  it  can  be  expected  to  exhibit  the  same  stability  and  economy.  In 


particular,  we  note  that  the  processing  times  for  this  algorithm  with  depth- 


dependent  background  velocity  will  be  conparable  to  those  for  a  constant 
background  k-f  algorithm.  In  addition,  we  shall  show  below  that  the 
algorithm  can  be  expected  to  be  quite  robust  even  when  the  "small” 
perturbation  assumption  is  violated. 

A  key  feature  of  our  approach  to  this  problem  is  repeated  application  of 
high  frequency  asymptotic  methods  to  obtain  an  inversion  formula  valid  in 
the  high  frequency  regime.  Discussion  of  the  motivation  and  justification 
for  such  high  frequency  approximation  may  be  found  in  BG  and  BCH.  In 
particular,  we  shall  use  a  ray  theoretic  Green's  function  in  formulating  our 
basic  integral  equation)  equation  (2)  below.  A  similar  approach  was 
presented  in  Clayton  and  Stolt  [1981].  Furthermore,  we  cannot  determine  an 
exact  inversion  of  the  fundamental  integral  equation  of  this  method. 
Instead,  we  assume  an  inversion  operator  which  consists  of  multiplication  of 
the  observed  data  by  a  factor  of  the  form  A  exp{-2iu>r}  and  integration  over 
the  data  set.  The  phase  t  is  the  traveltime  in  the  c(z)  background  medium 
between  the  source/receiver  point  and  the  output  point  at  depth.  Ve  must 
still  determine  the  amplitude.  A,  in  this  operator.  To  do  so,  we  require 
that  the  operator  applied  to  Kirchhoff  data  from  a  single  reflector  produce 
the  reflectivity  function  for  a  single  reflector,  to  leading  order, 
asymptotically.  This  leads  to  the  determination  of  the  amplitude  of  the 
inversion  integral  operator.  It  is  to  this  extent  that  we  then  claim  to 
have  an  inversion  operator  which  correctly  estimates  reflection  strength. 


-  2  - 


HIGH  FKEODENCT  STEP  LTD KAL  INVERSION 


Here  we  describe  the  formalism  for  determination  of  an  asymptotic 
inversion  operator,  np  to  the  amplitude.  A,  introduced  above.  We  employ  the 
same  wave  equation  model  as  described  in  detail  in  BCH.  If  v  is  the 
velocity  in  the  wave  equation,  we  set 


1 

v*(x,z) 


1 

c*(z) 


1  +  o(x,  z) 


x  *  (x,y)  . 


(1) 


Here  c(z)  is  the  known,  stratified  reference  velocity,  while  a(x»*)  i*  the 
desired  perturbation  correction  to  the  actual  velocity.  Furthermore,  we 
retain  the  assumption  of  backscatter  ("stacked”)  data.  In  this  case  the 
basic  integral  equation  for  a(x,z)  is  (cf.  BCH,  equation  (8)): 


a(x>z)  ^ 

—1 -  8  (x,*;Jjo») 

c  (z) 


£  -  U.n> 


(2) 


where  all  unmarked  integral  signs  are  over  (-<»,<») .  Here  u^  denotes  the 
backscattered  field  at  the  location  £  =  (£,n>  on  the  observation  plane,  i  = 

0  and  g  (the  "incident  field")  denotes  the  Green’s  function  corresponding  to 
the  stratification,  c(z).  In  contrast  to  the  constant  background  case,  g 
cannot  be  determined  exactly:  we  must  use  the  high  frequency  assumption. 

Fortunately,  this  assumption  is  completely  justified  on  the  geophysical 
exploration  scale  and  has  long  been  used  to  simplify  processing  formulas 
even  when  it  was  possible  to  derive  wide  band  analytic  results  (see  B(H) . 

We  use  J.  B.  Keller's  [1978]  ray  method  formalism  (see  also  Bleistein 
[1984]),  which  is  the  multi-dimensional  analogue  of  the  WKB  method  to  obtain 
a  parametric  representation  of  g  (see  Appendix  B) : 

. 


iwrCE, z) 


4n  \  k/E.O)  k}(E, z  )  E(E,z  )  H(E,z) 


Here,  if  we  introduce  the  transverse  distance. 


=  U  -  il  ' 


then  the  parameter,  E,  in  (3)  is  defined  as  a  function  of  p  and  z  by  the  ray 


equation: 


p  =  EHE,z) 


Further,  the  quantity  k^fE, z)  is  given  by 


(E, z)  =  Jn*(z)  -  E*  ,  (E*<  n* ( z) ) 


where  in  turn,  n,  the  index  of  refraction,  is 


n(z)  = 


"cTzT 


The  travel-time,  r,  is 


T  =  cTof  G(K 


,z)  .  G  =  . 

.0  k»(K’^ 


and  finally,  the  quantities  E  and  H  are  likewise  integrals  involving  n  and 
k( ,  These  integrals,  as  well  as  others  that  occur  subsequently,  are  defined 
in  Appendix  A.  There,  we  also  derive  some  needed  relations  involving  these 
quantities . 


Ve  shall  not  need  the  extension  of  k  in  (6)  to  the  range  n*(z)  >  K* 


because  onr  Fourier  transforms  are  integrals  over  real  wave  numbers  only. 


Thus  our  task  is  to  invert  the  integral  equation. 


us(£,<o) 


00 


a( x, z) 

c*  (z) 


2i<ot(K,z) 


k/K.O)  k,(K,z  )  E(K, z  )  H(K, z) 


(9) 


for  a(  x,  z)  in  terms  of  the  data,  u  (£.«)•  Again,  the  ray  parameter 
I  =  K(p,z)  is  defined  by  (5). 


Since  the  phase  in  (9)  resembles  that  of  a  Fourier  transform,  we  are 
motivated  to  seek  an  asymptotic  inversion  operator  of  the  form: 

W[f  (£,<■>) ]  (5' , x ' )  ~  //«•«/  e  2itan:(K'»z')  f(£,«)  ,  (10) 

with  p*  defined  as 

P'  =  |x*  -  £|  ,  (11) 

and  E'  defined  by 

p'  =  K'E(K' , z ’)  .  (12) 

Here  we  have  introduced  primes  to  avoid  confusion  with  the  integration 
variables  in  (9).  Applying  V  to  both  sides  of  (9)  and  writing  out  the  right 


hand  explicitly  we  have: 


WJUg] (**,**) 


-^T  J  dz  J  dw  |  |  d*x  J  J  d*£ 
16n  30  3  3  3  33 


exp{2iwi(K»  z, K' , z ' ) /c(0) } 

*  ,  .  “  A  P  ,z  k  ( K, 0)  k  (K,z)  E  (K,z)  H  (K.z>  * 

c  (z)  »  » 

where  we  have  introduced  the  phase , 

f  =  G(K,z)  -  G(K'.z')  . 

And  K,  K'  are  defined  by  (5)  and  (12)  respectively. 


(13) 


(14) 


Ve  perform  four  dimensional  stationary  phase  in  x  =  (x,y)  and  £  =  (£»t|) 
recognizing  that  K  depends  on  x  and  £  while  K'  depend.'  £.  Ht  present 
some  of  the  details  here  to  give  the  reader  the  flavor  or  the  type  of 
mathematical  analysis  involved.  Those  so  inclined  may  skip  directly  to  the 
result  given  by  equation  (21)  below.  Noting  that  (5)  and  (12)  allow 
computation  of  8K/dx,  3K/3y,  8K/dz,  etc.,  by  implicit  differentiation,  and 
using  the  result  of  Appendix  A,  we  find 


s  r  r 

§  =  G_  t —  =  G_ 

x  a  dx  K 


x— £ 

p(KE)K 


=  KB 


x-g 

pH 


x£ 

E 


(15) 


Similarly, 


! 


y 


*£  1 - 


where  we  have  introduced  the  short-hand 


E  =  E(K, z) ,  E’  =  E(K' , z)  . 


(16) 


(17) 


Thus  the  stationarity  conditions  are 


X  =  £  =  X ' 


(18) 


Bot  this  and  (5)  imply  (see  discussion  below  following  equation  48): 


K  =  0  . 


Simil iarly. 


K'  =  0 


These  results  greatly  simplify  the  remainder  of  the  stationary  phase 
calculation  (see  Appendix  D  for  the  details)  which  yields. 


W[ug]  (x ' . z ' )  ~  |  dz  a(x'.z)  A(O.z')  n(z)  j* 


1  d“  ®xp{2i7ToT  r.  n(t)dD  • 


It  should  be  realized  that  since  we  have  already  employed  the  assumption  of 
high  frequency  several  times  in  obtaining  (21)  and  shall  use  it  again  before 
we  finish,  the  u  integration  must  be  construed  as  a  band-limited 
integration.  These  matters  have  been  discussed  in  Cohen  and  Bleistein 
[1979b],  Bleistein  [1984],  and  BCH  [1984].  For  the  present,  we  shall  merely 
symbolize  the  effect  of  band-limiting  by  placing  a  subscript  B  on  the  Dirac 
delta  function  which  results  from  the  u-integrat ion  in  (21).  Thus, 


j"  d«o  exp[2i(<i>/c(0))  n(  df"}  =  nc(0)  n(t)df) 


=  nc(0) 


6b(z-z') 


and  so  (21)  becomes 
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A(0,i')o(x',i') 


(23) 


C 


a 


WIugHx',*')  ~ 


Solving  for  a,  writing  ont  the  expression  for  W[n  ](cf.  equation  10),  and 

0 


dropping  the  now  superfluous  prises  yields: 


a(x,z)  ~ 


nc(0) 


//*•«/ 


,  A(p,z)  -2iu>r(K,z)  , 

d“  TUTzT  e  ns(i.«o) 


with  K  determined  by 


p  =  KE  (K.z)  j  P  s  |*  ~  Si 


As  discussed  in  BCH ,  once  we  surrender  knowledge  of  the  low  frequency 
input  information,  we  cannot  obtain  output  trend  information.  It  is  to  be 
hoped  that  (by  iteration  if  necessary)  onr  c(z)  reference  velocity  is  an 
adequate  approximation  of  the  trend  to  the  depths  of  interest.  Vhat  we  can 
obtain  from  band-limited  information  is  a  perturbation  correction 
consistent  with  the  model  of  jumps  across  a  series  of  interfaces.  Ve 
determine  the  approximate  location  of  these  interfaces  as  well  as  the 
approximate  value  of  the  reflection  coefficient  at  the  interfaces.  This 
information  is  summed  up  in  the  reflectivity  function. 


P  =  ZRj6B(sj) 


where  s  is  a  (local)  arclength  variable  measured  normally  from  the  j 

J 

interface  and  R.  is  the  normal  reflection  coefficient  of  that  interface. 
J 

Clearly,  knowledge  of  £  is  equivalent  to  knowledge  of  reflector  location  and 
the  normal  reflection  coefficient  (see  equation  (43)  below).  In  turn  the 
latter  allows  direct  computation  of  the  jump  in  c  across  the  reflector. 
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According  to  the  theory  developed  in  Cohen  and  Bleistein  (1979b)  and 
reviewed  in  BCH ,  we  can  obtain  0  from  a  by  inserting  a  factor  of  iu/2c(z)  in 
(24)  to  obtain: 


~  nc<0?c(»y  I  i  B(p,l)  J  d“  “  e  2i"t(K'l)  (27) 

with  K  determined  by  (25).  Here,  since  the  inversion  depends  only  on  the 
ratio  of  A’s  in  (24),  we  have  introduced 


B( p , z) 


A(p,z) 

TUT,  IT  * 


(28) 


However,  and  this  is  the  key  point:  at  this  stage  of  the  derivation  we  have 
no  information  about  B(p,z)l  Any  choice  of  this  quantity  gives,  in  the 
language  of  BCH,  a  structural  inversion,  i.e.,  a  migration.  In  order  to 
determine  a  choice  of  B  which  will  yield  an  accurate  approximation  of  the 
interface  reflection  coefficients,  we  will  insert  in  (27),  a  canonical  set 
of  scattering  data,  u  ,  and  then  determine  B  by  enforcing  the  asymptotic 
equality  in  (27).  For  this  purpose  we  will  use  Kirchhoff  data  for  a  single 


reflecting  surface. 


a 

1 


« 


€ 


p 


i 


4 

r 


DETERMINATION  OF  TD  INVERSION  AMPLITUDE 

In  order  to  find  an  amplitude  B(p,z)  for  the  inversion  f omuls  (27),  ve 
first  obtain  an  expression  for  the  Kirchhoff  representation  of  data  Ug. 
Such  data  employs  the  high  frequency  assumption  (by  using  the  multi¬ 
dimensional  WKB  representation  of  the  incident,  reflected  and  transmitted 
fields),  but  does  not  make  the  Born  approximation  of  small  reflection 
coefficient.  Thus,  if  ve  can  determine  a  B  vhich  enforces  asymptotic 
equality  in  (27)  for  such  data,  our  algorithm  is  likely  to  be  quite  robust 
for  large  contrast  interfaces. 

It  remains  to  decide  on  the  surfaces  to  use  in  computing  the  Kirchhoff 
approximation  to  u  .  It  turns  out  that  a  single  surface,  such  as  the  tilted 
plane,  z  =  zQ  -  x  tan  0,  vould  suffice  to  determine  B.  However,  we  carry 
out  our  calculations  for  the  more  general  cylindrical  (i.e.  y  independent) 
surface : 

z  «  h(x)  .  (29) 

One  may  think  of  the  tilt-plane  members  (which  are  included  in  (29))  as 
determining  B,  while  the  remaining  members  conf irm  this  choice  of  B. 


In  Appendix  C,  we  show  that  the  back  scatter  at  £  from  (29)  has  the 
Kirchhoff  (high  frequency)  representation: 


Ug(£,u) 


~  2iu> 


y  RS 


2iwt(K,l) 

e 


J 


(30) 


subject  to 


4 
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P  «  EE  (K,i)  ,  p  *  |l  -  jj,  z  *  h(x) 


Here  we  here  need  bars  to  distinguish  the  spatial  variable  from  the  output 
variables  in  the  inversion  formula  (27)  and  have  also  introduced  the 
quantities : 


f  J2=S2±1  -  k  (ifI)  1 
L  E(f,«)  1  i 

[  c(0)  l|l+h'*  ] 


16n  k,(K,0)kt(E,Z)E(U)l(K,Z) 


and  the  (non-normal)  reflection  coefficient* 


’  Yx“*8n<T)  h+-r-^r  • 

J  Cj  c 


In  turn,  ct(z)  denotes  the  actual  velocity  below  the  reflector,  z  =  h(l) , 


that  is. 


c^z)  *  c(z)  4  Ac 


where  Ac  is  not  accounted  for  by  the  stratified  reference  profile  c(z). 
Obviously  determining  R  is  tantamount  to  determining  Ac  or  c1(z).  Combining 
(27)  and  (28),  we  have 


Ki-*’  ~  ,.(o‘c(tS  J  <■»«’/  J  i‘?  f  f 

rRSB  e21“  [  ] 


Our  goal  is  to  determine  B  so  as  to  obtain  the  asymptotic  equality. 


We  now  carry  out  stationary  phase  in  I,  f,  q  (see  Appendix  D  for 
details)  and  obtain  the  stationarity  conditions: 


These  imply: 


y  =  q  =  y 


X  -  5  =  -sgn(h’)K  E(K.Z) 


x  -  5  *  -  sgn(h')  KE(K.z) 


i .  i  .  fliUiiL  . 


^(K.I)  = 


where  again  X  -  h(x) . 


Geometrically,  these  conditions  confirm  the  cylindrical  nature  of  our 
reflector  and  show  that  the  output  point  (x.z)  lies  on  a  specular  ray. 
Furthermore  (37-39)  yield 
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(42) 


r  "  "7 


*x  =  -T 


which,  in  turn,  imply  that  R  reduces  to  the  normal  reflection  coefficient. 


ci  -  c 

R  -  Rn  -  7— TT 


Completing  the  stationary  phase  analysis  (again  see  Appendix  D)  we  find 


16nn  (z)  i  1+h*  RnSB(p,z)  f  2i«[T(f.l)  -  t(f.z)] 

p(x.z) - — . - . — -  ‘  due  (44) 

c(0)  J Idetl, . i  J 


Here, 


1  (l-fh,a)*  t  n^h*  T  1 _ 1  1 

(K.z)E(K.Z)  H(K.z)H(K,Z)  ^  +  k,s  [  H(f,z)  H(K,Z)  J 


However,  the  final  integration  in  u,  yields  a  delta  function  whose  argument 


can  be  transformed  to  arclength  along  rays: 


J  dm  e2i«[r(K,Z)  -  t(I.x)]  =  nbglrd.Z)  -  r(t,z)) 


=  nc(0)  6bIG(I,z)  -  G(K.z)] 


^(K.z) 

=  nc(0)  — -  SB[z  ~l] 

n*(z)  B 


Jtc(O)  •  6g[s(z)  -  s(z)l 


Here  the  last  equality  involving  the  ray  arclength  variable  follows  from 
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8)  and  Appendix  B.  Ve  say  note  that  since  z  «=  I  when  the  delta  fraction 


"acts."  the  stationarity  condition  (38-39)  imply  that  the  output  point  (x,z) 
coincides  with  the  specular  point  ( x ,  z)  =  (x,  h(D).  Furthermore,  when 
z  -  I  the  second  term  in  the  square  brackets  in  (45)  drops  out  and  hence 
using  the  definition  (33).  we  see  that 


16jt*(l+h*  *)  kt  (1,0)  kf  (K.z) 


(47) 


Hence  (44)  reduces  to 

R  n(Z)  B 

p(x,z) . . . — - - 6[s(z)  -  s(z)] 

)|l+h'*  k,  (1,0)  k/t.z) 

R  B 

’■E7?or“*(',-*<,),  <4,) 

R  B 

=  -  n  8  8[s(z)  -  s ( z ) ] 

)|l-E* 

Here,  the  middle  equality  follows  from  (41),  while  the  final  inequality 
follows  from  the  definitions  (6-7)  when  z  =  0.  Also  the  quantity  B#  denotes 
the  value  of  B  at  the  stationary  point  and  when  z  =  h(x).  Thus  to  enforce 
asymptotic  equality  in  our  inversion  algorithm  (27),  we  need  to  choose  a  B 
which  reduces  to  ljl-1*  under  these  conditions:  Obviously  such  a  B  is 


B(p,z) 


(49) 


where,  as  usual,  the  parameter  E  is  determined  by  (25).  Even  though  the 
algorithm  (27)  depends  only  on  B,  and  not  on  the  A(p,z)  originally 


introduced  in  (10),  self-consistency  demands  that  ve  be  able  to  write  B  in 
the  fora  (28).  This  is  indeed  simple  to  do  since,  for  example,  if  we  pick 


A(p,z) 

Then  by  (25),  and  (A-7)  and  (A-15) 
K  =*  0.  Thus, 


the  unique  E  corresponding  to  p 


(50) 

=  0  is 


A(0,z)  *  1 


(51) 


and  (49)  verifies  (28).  The  fact  that  A  is  non-unique  is  irrelevant  since 
(27)  only  depends  on  the  ratio  B. 


It  should  be  noted  that  choosing  B  as  in  (49)  only  guarantees  that  (27) 
gives  an  accurate  estimate  of  0  at  the  specular  point  (stationary  point). 
While  this  is  asymptotically  the  most  important  point,  any  practical 
implementation  of  (27)  involves  integrating  over  a  region  including  the 
specular  point.  It  is  conceivable  that  a  still  more  accurate  B  could  be 
discovered  by  somehow  determining  it  without  use  of  the  canonical  Kirchhoff 
data.  We  have  made  a  considerable  effort  to  find  a  suitable  generalization 
of  Fourier  inversion  which  would  allow  a  derivation  along  the  lines 
originally  presented  in  Cohen  and  Bleistein  (1979a)  and  reviewed  in  BCH . 
Unfortunately,  to  date,  these  efforts  have  not  succeeded. 


The  algorithm  for  the  reflectivity  function  derived  in  the  previous 
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i 


s 

* 


: 


sections  is 


P(x,z)  ~ 


8i 


nc(0)  c 


^-jjA 


-2iwc(K,  z) 


•  dt  Us(t.$)  • 
^0 


iut 


(52) 


p  =  KE(K,z)  ,  P  s  |x  -  £| 


Here  the  integrals  E  and  G  are  defined  in  Appendix  A  and  Ug  is  the 
backscatter  data  observed  on  an  areal  array.  For  actual  data  processing,  it 
is  convenient  to  "fold”  the  unphysical  negative  frequencies  onto  the 
positive  ones  by  replacing  to  by  -»  on  the  interval  (-*,0).  At  the  same 
time,  we  introduce  the  physical  frequency  variable  (measured  in  Hz): 


(53) 


and  explicitly  acknowledge  the  bandlimiting  by  introducing  F(f),  a  tapered 
high  pass  filter.  After  these  changes,  we  have: 


&(x,z) 


~64n 

cToTcTzT 


•  In  f  df  f  F(f)  e~4rrifT(K,z) 
J0 


p  =  IE  (k,z) 


f  dt  Us(t,£)  e2vift 

J0  * 

P  ■  Is  -  £l  • 


(54) 


__  1  C  __ 


In  practice,  areal  observations  are  often  not  available  and  instead  only 
a  linear  set  of  data  is  nsed.  In  this  case  ve  cannot  hope  to  reconstruct  a 
three  dimensional  image  of  the  snbsnrface  and  instead  seek  a  two  dimensional 
slice,  £(x,0,z)  ■  0(x,z),  consistent  with  the  data  available.  Since  the 
data  is  now  independent  of  n, 

Us<t.£)  =  Us(t,$.0)  ■  Ds(t,?)  ,  (55) 

and  we  may  carry  ont  an  additional  stationary  phase  calculation  in  i).  The 
stationarity  condition  is 


H  =  y  (56) 

and  the  analogue  of  (52)  is  found  to  be  (see  Appendix  D  for  details): 


dt  Us(t,?)  eiWt 

■*0 

p  *  KE(I,z)  .  p  -  |x-C  |  . 


(57) 


while  (54)  becomes: 
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p(x,z)  ~  -  32n- -  d5  J(l-K*)  E(K, z)  . 

^ToT  c ( z )  J 

(Re  -  Im)  f  df  J?F(f) 

J0  1 

(58) 

•  f  dt  U_(t,$)  e2"ift  . 

J0  b 

p  =  KE  (K.z)  .  p  =  |x-5  | 

The  basic  concepts  of  reducing  (58)  to  a  computer  code  are  the  same  as 
those  discussed  in  BG  for  the  algorithm  presented  there.  Briefly,  the  t  and 
f  integrals  are  performed  routinely  using  an  efficient  FFT  algorithm.  The 
main  complication  in  (58)  lies  in  the  expressions  E(K,z)  and  t(K,z),  both 
being  integrals  defined  by  (A-4)  and  (8)  respectively.  This  is  a  bit  subtle 
in  that  the  parameter  K  (see  Appendix  B)  can  be  viewed  as  determining  the 
starting  angle  for  a  ray  connecting  the  surface  point  ((,0)  to  data  point 
(x,z)  in  (58).  Therefore,  for  a  given  offset  p  =*  |x-{|,  K  is  defined  by  the 
implicit  relation  p  *  KE(K,z).  In  the  computation  this  issue  is  handled 
quite  efficiently  by  two  tables  for  evaluating  t(K,z)  and  the  amplitude 
(involving  E)  in  (58)  as  functions  of  p  and  z. 

The  computation  time  of  the  resulting  algorithm,  as  pointed  out  in  BG,  is 
comparable  to  a  standard  k-f  migration  algorithm  for  constant  coefficients. 

Ve  are  now  developing  and  testing  a  very  carefully  designed  Fortran  77 
computer  code  for  (58)  and  expect  to  present  the  results  of  those  tests  in 


the  near  future. 


CONCLUSIONS 


Ve  have  presented  the  derivation  of  an  inversion  algorithm  for 
backscattered  ('stacked')  seismic  data.  We  made  four  major  assumptions: 
(i)  the  acoustic  wave  equation  is  an  adequate  model,  (ii)  backscattered 
data  has  amplitude  information  vorth  preserving  fairly  accurately,  (iii)  the 
actual  reflectivity  coefficients  can  be  adequately  modeled  as  perturbations 
from  a  continuous  reference  velocity  which  depends  only  on  the  depth 
variable,  (iv)  the  subsurface  can  be  adequately  modeled  as  a  series  of 
layers  with  jump  discontinuities  in  the  velocity  (or  impedance)  at  these 
layers. 

The  last  assumption  is  unavoidable  given  the  nature  of  the  high  pass  (on 
the  exploration  scale)  data  collected  in  the  field.  The  third  assumption  is 
inherent  in  our  approach  although,  as  pointed  out  above,  the  algorithm  can 
be  expected  to  be  robust  even  when  this  assumption  is  violated.  Also  the 
algorithm  presented  here  represents  a  considerable  improvement  over  earlier 
algorithms,  such  as  Cohen  and  Bleistein  [1979a] ,  which  perturbed  from  a 
constant  reference  velocity. 

On  the  other  hand,  weakening  of  the  first  two  assumptions  seems  eminently 
feasible  and  we  hope  to  apply  the  techniques  expanded  in  this  article  to 
both  inversion  of  offset  data  ('inversion  before  stack*)  and  to  equations 
which  more  accurately  describe  the  wave  propagation  in  the  earth. 

It  is  already  known  (see  BG)  that  algorithms  with  the  structure  of  the 
one  presented  here  are  numerically  stable  and  are  computationally  efficient 


relative  to  other  seismic  data  processing  algorithms. 


The  essence  of  the  derivation  presented  is  to  form  an  inversion  operator 
with  a  suitable  pre-specif ied  phase,  bnt  unknown  amplitude.  The  amplitude 
is  then  determined  by  demanding  that  for  high  frequency  synthetic  scattering 
data  from  a  fairly  general  surface,  the  algorithm  asymptotically  produce  the 
correct  reflectivity. 
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APPENDIX  A 


We  define 


and  the  integrals. 


NOTATIONS  AND  IDENTITIES 


n(z) 


c(0) 
c  ( z )  ' 


(A-l) 


(K, z) 


'|n*(z)  -  K*  , 


(A-2) 


D(K.z) 


(A-3) 


E(K, z) 


(A-4) 


F(K.z) 


(A-5) 


G(K,z) 


(A-6) 
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H(E,z) 


(A-7) 


n*(Ddf 

^(E.f) 


s(K,z) 


(A-8) 


Similar  quantities  occur  in  the  T-p  theory,  see  Diebold  and  Stoffa 
[1981].  Among  the  many  relations  which  link  these  quantities,  we  cite  below 
those  that  are  useful  in  carrying  out  the  calculations  presented  in  this 
paper  and  its  appendices. 


First  of  all,  from  (A-2)  it  follows  that 

D  +  K*E  =  G  , 


(A-9) 


E  +  E*F  =  H  .  (A-10) 

Next  we  cite  the  k  and  z  partial  derivatives  of  D,  E  and  G  which  follow 
respectively  from  use  of 


and  the  Fundamental  Theorem  of  calculus: 

D_  =  -  KE  .  D  -  k  (E,  z)  , 
K  z  > 


( A-ll) 


(A-12) 
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E„  =  KF  .  E  =  -A- 

z  *» 


(A-13) 


Gr  =  KH  ,  G  =  ■£— 
K  z  k 


(A-14) 


Finally,  from  (A-13)  and  (A-9)  it  follows  that 


(KE)k  =  H  . 


(A-15) 
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APPENDIX  B 


THE  S1BAXIFIED  MEDIA  GREEN'S  FUNCTION 

Using  Keller's  [1978]  "ray  method*,  developed  in  the  1950's,  ve  seek  s 
high  frequency  spprozimation, 

g(z.z)  ~  A(z,z)  (  j  =  (z.y)  (B-l) 

which  ssymptoticslly  satisfies  the  Helmholtz  equation. 

» 

V*g  +  — r- —  8  =  -8(z-J)  6(z)  ,  £  =  (?,i\)  .  (B-2) 

c  (z) 

To  complete  the  specification  of  g.  we  insist  that  it  behave  like  the  free 
space  (i.e.  constant  c)  Green's  function  as  the  field  point,  (z.z), 
approaches  the  source  point,  (£.0) .  This  entails  the  conditions, 

t  -4R/c(0)  .  A  -»~g-  (B-3) 


as  R  -4  0,  where 


R*  =  |z-£|*  +  z*  .  (B-4) 

Ve  substitute  (B-l)  into  (B-2)  and  separately  equate  the  coefficients  of 
and  u  to  zero  (this  is  the  high  frequency  spprozimation)  giving  rise  to 
the  eikonal  equation, 

k  •  k  =  n  ( z)  ,  k  =  c(0)Vr  ,  n  =  ( B-5 ) 


and  the  transport  equation. 


2  k  •  VA  +  (V  •  k)  A  «  0  . 


(B-A) 


The  former  equation  can  be  solved  by  the  method  of  characteristics  (see 
Bleistein,  1984)  which  reduces  the  problem  to  the  eolation  of  a  system  of 
ordinary  differential  equations.  The  first  of  these  equations  are: 

dz  dz 

^  =  l  •  SSm  k.  *  IB<kx*ks>  »  l  -  <I*k,>  <B-7) 


d|  dk, 

•j—  *  0  .  -rrr~  m  ~  n  n*  (z)  (B-8) 


which  define  the  rays»  o  being  the  ray  parameter.  The  source  term  in  (B-2) 
makes 


J(0)  *  5  ,  z(0)  =  0  .  (B-9) 

a  natural  choice  as  initial  data  for  (B-7) .  The  data  for  k(0)  consists  of 
an  arbitrary  unit  vector  (cf.  (B-5) ,  noting  that  n(0)  =1).  To  be  specific 
we  choose 


5(0)  arbitrary  ,  k,(0)  = 

where  we  have  introduced 


(B-10) 


K  =  III  =  \  kI  +  kI  • 


(B-ll) 


From  (B-10)  and  the  fact  that  we  are  not  considering  turned  rays  here,  it 
follows  that  K  =  (k4,  k2)  can  be  viewed  as  the  direction  numbers  of  the  rays 
initiating  from  ({.  q,  0). 


■*"<r  ■  r — r — *  ■ .  w  ;  ■ 


^  ^  ■■  i  - 


p 


■-  ■'V 


Using  (B-6)  snd  (B-7)  we  find 


dt  n 
do  *  c(0) 


2  “  *  (V  •  k)  A  =  0 


(B-12) 


(B-13) 


which  together  with  the  date  (B-3)  conpletes  the  specification  of  g, 
expressed  in  (B-l),  in  terns  of  a  system  of  ordinary  differential  equations. 

Proceeding  to  the  analysis,  we  first  note  that  by  (B-8),  K  =  K(0),  so 
henceforth  we  simply  write  K  for  this  constant  vector.  Then  the  eikonal 
equation  (B— 5)  gives  ns 


c.U.z)  =  Jn*  -  i  . 


(B-14) 


Then  (B-7)  yields 


dx 

~dz 


'In*  -  K1 


(B-15) 


or 


i  -  i  »  K  E(I,z) 


where  E  is  defined  by  (A-3) .  Similarly,  we  find 


(B-16) 


1 

tutt 


X  *  -r^T-  O(K.x) 


(B— 17 ) 


where  0  is  defined  by  (A-6). 


•.  ■ - 
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If  we  introduce  the  ray  Jacobian 


d(x.z) 

3(K,o) 


(B-18) 


and  use  the  fact  that 


=  J  V  •  k  , 


(B-19) 


then  the  transport  equation  (B-6)  can  be  recast  as 


constant 


(B-20) 


A  Jj  =  lim  A  l|7 

R->  0 


(B-21) 


which  by  (B-3)  implies  that 


J7R-.0  [<’* 


(B-22) 


Ve  now  indicate  how  to  obtain  the  partial  derivatives  which  are  the 
elements  of  J. 


Implicit  differentiation  of  (B-17)  with  respect  to  ka>  ka  and  (A-14) 
yields 
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a  k  M 

—  -  *  j  1=12 
0k  a  *  J  1,2 

J  D 


(B-23) 


Similarly  from  (B-16) 


Zi.Et  .  [p.A]  k  k  .  j  .1.2 

j  n  J 


(B-24) 


dz  3z 

Since  and  are  given  directly  by  (B-7),  we  are  now  able  to  form  J. 


A  short  calculation,  involving  the  ose  of  (A-10)  yields 


J  =  k  EH  . 
a 


It  is  easy  to  show  that  as  R  -4 0  (equivalently  a  -40) 


In  the  same  limit.  (B-16)  implies  that 


(B-25) 


(B-26) 


kj-4  kj  (K.O)  =  <1-K  =  z/R 


(B-27) 


so  that 


TTOT  • 


(B-28) 


Then  (B-22)  gives 
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1 


(B-29) 


A  =  —  -  •  - '  . 

Jks(K,z)E(K.s)H(E.x)  4n  jFTOT 

Since  A  end  t  depend  only  on  K,  we  need  only  employ  the  magnitude  of  (B-16) 
in  the  sequel : 

p  =  KE(E,z)  j  p  *  |x  -  $|  .  (B-30) 

Equation  (B-14) ,  (B-17),  (B-29)  and  (B-30)  are  equivalent  to  equations  (3-8) 
of  the  text. 

Finally,  we  relate  our  parameter  a  along  the  rays  to  the  corresponding 
arclength  parameter.  From  (B-7)  we  have 

dx  dx  dz  dz 

-  .  -  + - =  K  +ks=n  (B-31) 

dcr  do  do  do 


and  so 


IT-  »<*> 

do 


Using  the  second  equation  in  (B-7)  once  more,  we  find 


ds _ n_ 

dz  l. 


or 


n(  d 

k^K.t) 

0 


(B-32) 


(B-33) 


(B-34) 
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APPENDIX  C 


mCHHOFF  DATA 


The  derivation  presented  in  Cohen  and  Bleistein  [1983]  applies  here  with 
the  constant  c  Green's  fnnction  used  there  being  replaced  with  the  c(z) 
Green's  function  derived  in  Appendix  B.  Thus,  equation  (8)  of  that  paper 
stay  be  recast  in  our  present  notation  as 


Ug(z,o>)  ~  |  dS  R-—  g*(x,z;J 


(C-l) 


with  g  defined  by  equations  (3-6)  and  the  reflection  coefficient  R  being 
defined  by 


r  -  rx 
r  +  r± 


(C-2) 


where 


r  =  #  •  Vr  ,  rx  =  sgn(y)  I  r*  +  —  -  •  * 


c,  c 


(C-3) 


and  8  ia  the  upward  normal  to  S.  Here  ca  is  defined  by  (35). 


Since  g  has  the  form  given  in  (B-l) 


g*  ~  2iw  8  •  Vt  g*  , 


(C-4) 


subject  to 
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p  =  K  E (K, z) 


( C-5) 


Bence  (B-l),  with  the  definition  of  y  given  in  (C-3)  and  that  of  S  given  by 
(33),  yields 


3  a 

7?  8 


2iu>yS 


„2ill)T 

e 


(C-6) 


Thus  (30)  follows  from  the  form,  (29),  of  the  surface.  It  remains  to 
establish  the  detailed  calculation  of  y  given  in  (32). 


We  have: 


y  =  8  •  Vt 

=  TToT  8  *  VG(K’z) 


( C-7) 


1  (h ' ,0,-1)  r  „  3K  „  3K  „  dK  .  „  1 

cW  V-  -  -  ~  l  G*  37  '  gk  W  K  37  +  Gz  J 


{l+(h‘ 


) 


1  1  r  3K  „  3K  „  1 

c(0)  , - r  l  h  GK  3x  GK  3z  ~  Gz  J 

)|l+(h') 


From  the  constraint  (C-5)  we  compute 


3K  x-5 
3x  “  p(KE)K 


3K 


EE 

z 


3K  =  "  IK  E) 


( C-8) 


and  then  (32)  follows  from  (C-7),  (C-8)  and  the  results  of  Appendix  A. 


APPENDIX  D 


STATIONARY  PEASE  CALCULATIONS 

Assuming  that  the  phase,  |(z),  has  a  single  simple  stationary  point,  z#: 

Vf(z  )  =  0  .  det  i  (z  )  *0  (D-l) 

-s  -  zizj  s 

the  integral, 

I(X)  ~  J  eU,(l)  A(z)d”z  ,  X>0  (D-2) 

can  be  evaluated  asymptotically  by  the  multidimensional  stationary  phase 
formula,  (see  Bleistein  [1984]  or  Bleistein  and  Handelsman  [1975]). 

I(X)  ~  [4^-1  - - - ezpfiXl  +i-4*ig*,J  »  (D-3) 

1  J  IP^TT 

Here  A#>  |$,  and  ijj*  denote  respectively  the  amplitude  A,  the  phase  i 
and  the  Hessian  matriz,  fx  x  evaluated  at  the  stationary  point,  z  =  z#. 
Further,  sig  fjj  denotes  the  signature  of  l^j,  i.e.  the  number  of  positive 
eigenvalues  minus  the  number  of  negative  eigenvalues. 

If  there  were  several  stationary  points,  (D-5)  would  be  replaced  by  a  sum 
over  the  contributions  from  these  several  points.  If  the  Hessian  matriz 
vanished,  (D-3)  would  have  to  be  replaced  by  a  more  general  result. 
However,  only  (D-3)  is  required  here. 

In  our  applications,  instead  of  a  positive  parameter,  X,  we  have  the 
signed  quantity,  2w/c(0).  The  sign  can  be  accounted  for  by  replacing  i 
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by  -  i.  la  addition,  only  the  caae  n  =  4  occur*  in  the  text,  hence 


i*  A 


>'*>  -  [tit1] 


dett . . 
ij 


+  -p  sgn(w)  sig  fj  }  .  (D-4) 


The  context  for  our  asymptotic  evaluations  is 

J  =  |  I(«)  u  dfa> 

(more  precisely  a  handlimited  version  of  (D-5)) 


(D-5) 


The  most  important  case  for  us  is  when  sigi^j  =  0,  so  that 


r  i*  * 

1  nc(0)  I  - 5 - 

■  .  [  ».  ] 

L  J  J  detlT 

[  ig  ]  .  ( sigilj  =  0)  ,  (D-6) 


where  6  denotes  the  Dirac  delta  function.  The  only  other  case  that  occurs 
below  is  sigfjj,  =  ±2,  in  which  case, 

*  2  I111  f 

IU)  ~  [  ~T~r"  ]  *  *  *gn(<o))  e  *  ,  (sigijj  =  -2)  ,  (D-7) 


n*c* (0) 


1 

T  ’ 


(*ig*ij  =  ±2)  * 


(D-8) 


Although  the  evaluation  of  the  signature  of  a  symbolic  4x4  matrix  can 
be  tedious  or  impossible,  our  task  is  considerably  simplified  by  the  fact 


that  our  Hessian  matrices  have  the  special  form. 


« 


T 

whence 

r« 


« = 


a 

0 


v 

0 


0  v  0 
0  0  u 
0  y  0 
u  0  6 


det  i  =  (ay  - 


v  ) 


(06  -  u  ) 


(D-9) 


( D-10) 


€ 


To  evaluate  the  signature,  we  need  to  determine  the  roots,  X,  of  the 
eigenvalue  equation. 


det  (i  -  XI)  =  0  » 


( D-ll) 


where  I  is  the  identity  matrix.  From  (D-10),  we  find  at  once: 


det  (|  -  XI)  =  [  X2  -  ( a+y)  X  +  ay  -  v*  J  [  X*  -  (0+5)  +  06  -  u*  J  .  (D- 


12) 


Thus , 


2  2 

v  >  ay  ,  u  >  06  =>  sigl^  =  0 


( D-13 ) 


and  similarly. 


v  >  ay  ,  u  <  05  =>  sigf^  =  • 


(D-14) 


etc. 


We  now  turn  to  the  specific  stationary  phase  calculations  in  the  text. 
We  shall  freely  use  the  results  of  Appendix  A  without  explicit  citation. 
First  we  examine  the  phase  in  (14): 
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§<x,£)  =  G(K, z)  -  G(K'.z') 


where  I(p,i),  K'(p',z')  are  defined  by 


p  =  KE  (K,  z)  .  p'  =  K'E(K'.z') 


with 


p  -  I*  -  il  •  p'  =  Is'  "  £l  • 


Implicit  differentiation  of  (D-16)  with  respect  to  x^  yields 


x.  -  $.  dK 

i.  =  (KE)K-5— = 

i  i 


so 


aK 


xi -  5i 


3x.  ~  KE(K,z)  IKK, z)  ' 


Similarly, 


dK 


h  -  Xi 


3K’ 


h  -  *i 


Yta.zwu.z)  •  -  nraira 


Now  using  these  resalts,  we  find  that 


>  _  „  3K  _  Xi  '  5i 

#xi  -  gk  al^  -  TTTTzT 


(D-15) 


and  similarly 


*?i  ~  E(K.z)  E( K' , z ' )  * 

From  (D-21)  and  (D-22),  we  see  that  the  stationarity  conditions  are: 

x  =  £  =  x’  . 

(D-23 

Hence  by  (D-15) > 

—  2  — 

p  =  p'  =  0  . 

(D-24 

Now  it  is  clear 

that 

K  =  K'  =0 

(D-25 

is  a  solution  of 

the  constraint  conditions  (D-16) . 

Moreover,  since 

(ke>k  =  H  >  0  . 

(D-26 

this  is  the  only 

solution.  We  obtain  for  the  phase 

at  the  stationary  point 

i  =  G(0.z)  -  G(O.z’)  -  f  n(t)dt 
s  J  . 

(D-27 

z 

and  for  the  components  of  the  Hessian  matrix: 

1 


a  -  1/E(0,z)  ,  o'  -  1/E(0,z')  . 

In  the  notation  of  (D-7)  we  have 

a  a 

v  -  oy  =  n  -  p8  =  oo'  >  0 


(D-29) 


D-30) 


so  that 


_  _a 

d6t  *ij  =  [  E(0,z)  E(O.z')  ]  '  ,ig#ij  =  ° 

These  results  allow  ns  to  nse  (D-4)  in  (13)  to  obtain 
obtain  ( 22) . 


(D-31) 

(21)  and  (D-6)  to 


We  now  turn  to  the  stationary  phase  evaluation  of  (36).  Here  the  phase 
is 


i  «  G(K,X)  -  G(K, z)  ,  (D-32) 

subject  to  the  usual  constraints 

P  -  I*  ”  4l  “  KE(K.z)  .  J5  =  |x  -  £|  =  1  E(f ,X)  .  (D-33) 

This  stationary  phase  evaluation  is  somewhat  more  difficult  because  X 
depends  on  X: 

X  *  h(x)  ,  (D-34) 

and  hence  the  |j  condition  is  more  complicated  than  that  in  (D-21) .  The 
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analogue  of  (D-18)  is 


-  (KE)_  —  +  KE  »  H  —  +  ~  h'  .  (D-35) 

P  dt  Z  3X 


at  s  x-£  _  fv_ 

dx  KEH  kaH 


where  we  use  notations  lit) 


1  «  E(i,X)  (D-37) 

for  quantities  which  depend  on  1  and  X.  The  remaining  partials  of  K  and  K 
have  the  same  form  as  in  (D-18)  and  (D-19)  since  X  does  not  depend  on  ^  Jj, 
or  y,  and  since  a  is  a  completely  independent  variable: 

il  =  23  .  il.ii  .  ...  (d-38) 

ay  KEH  dt  KEH 


These  results  allow  one  to  compute: 

«.  h,<,) 

=  izi  _  ll  h- 


!  _-23.  .  f  ,  f  =  £11  -  JTX 

7  E  ?EE  nEE 


(D-3 9) 


The  y»  T)  derivatives  give  rise  to  the  single  stationarity  conditions: 
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y  -  n  *=  y 


(D-40) 


Hence  the  constraints  (D-30)  reduce  to 


P  -  |x-?|  =  EE(I,z)  .  p  =  |x~?|  =  lE(I.x)  , 


(H-41) 


which  can  be  rewritten  as 


x-{  =  pKE  ,  p  ■  sgn(x-$) 


(D-42) 


=  pSE  ,  p  ®  sgn(x~5) 


(D-43 ) 


These  results  allow  us  to  simplify  the  renaining  partials  in  (D-39)  to 


fx  *  pi  +  k4h' 


(D-44) 


#y  =  -  pi  +  pi  . 


(D-45) 


For  stationarity.  we  nust  have 


p  «  p  *  -  sgn(h') 


and  then  also 


(D-46) 


1  =  1  =  kj  h '  . 


(H-47) 


The  latter  equality  allows  to  deternine  that 
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(D-48) 


i  ,  »<1)  |h'j  #  . 

Jl+h'*  Jl+h'* 

Finally  (D-46)  allows  ns  to  restate  (D-42),  (D-43)  as 

x  -  5  =  -  1  E(I, z)  sgn(h') 

and 

x  -  5  =  -KE(i,5)  sgn(h') 

At  the  stationary  point,  the  phase  is 

fs  -  G(l.l)  -  G(f,z)  . 


Further,  a  number  of  terms  of  the  Hessian  matrix  vanish  because  of 
leaving  us  with 


0 

*t5 

0 

0 

i__ 

yy 

0 

*yn 

•h 

0 

0 

0 

** 

0 

which  has  the  form  (D-9)  with 


t 

v  -  ay 


-  *«  »«>l. 


and 


-  Pr  =  <«'  -!-_#) 

yn  yy  nn 


(D-49) 

(D-50) 

(D-51) 

(D-40) 

( I>— 52) 

(D-53) 

(D-54) 


Ve  observe  that  (D-36)  and  the  stationarity  conditions  yield 


||.|  ‘llilL.idrt’*)  .  n--„»U') 

"  Is  B  B  B 

and  similarly 

1£  I  -  £  H  I  -P 

»«  Is  "  S  '  9«  Is  "  e<i,«)  ' 

These  facts  allow  to  evaluate  (D-53)  as 

*  (1+h'V  .  n(Z)h"  r  1  11 

v  -  =  — - = —  +  7 - —I - rr 

H(K, z)H(K,  z)  J1+h,s  Lh(K.z)  H(K,z)J 

and  (D-54)  as 

u*  -  py  - - - - >  o  . 

E(I,z)E(K,z) 


(D-55) 


(D-56) 


(D-57) 


(D-58) 


Thus,  "near"  the  reflector,  z  =  I.  both  v*>aY»  u*>p6  hold  and  we  have 
sigi^j  =  0.  In  this  regime,  (D-6)  applies  and  we  see  that  we  have  a  delta 
function  "acting”  at  the  reflector.  On  the  other  hand,  far  from  the 
reflector,  we  may  have  v*<uy  leading  to  the  distribution,  (D-8) .  However, 
since  z  is  "far”  from  z,  this  distribution  is  regular,  indeed  "small”,  so 
that  the  results  given  in  the  text  follow. 
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